INVERSE SCATTERING IN MULTIMODE STRUCTURES 
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Abstract. Wc consider the inverse scattering problem associated with any number of interact- 
ing modes in one-dimensional structures. The coupling between the modes is contradirectional in 
addition to codirectional, and may be distributed continuously or in discrete points. The local cou- 
pling as a function of position is obtained from reflection data using a layer-stripping type method, 
and the separate identification of the contradirectional and codirectional coupling is obtained using 
matrix factorization. Ambiguities are discussed in detail, and different a priori information that can 
resolve the ambiguities is suggested. The method is exemplified by applications to multimode optical 
waveguides with quasi-periodical perturbations. 
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1. Introduction. In waveguides that support several modes, scattering, or cou- 
pling between the different modes, may appear due to different kinds of perturbations. 
Possible perturbations are reflectors, gratings, bends, tapering, and other kind of ge- 
ometrical or material modulation along the waveguide. The coupling may be both 
codirectional (coupling between modes that propagate in the same direction) or con- 
tradirectional (coupling between modes that propagate in opposite directions). The 
direct scattering problem of computing the scattered field when the probing waves 
and the scattering structure are known has been extensively discussed in the liter- 
ature [211 [38l [21]. The inverse scattering problem associated with two interacting 
modes is also well understood, and has been treated in several contexts since the 
C\| [ pioneering work by Gel'fand and Levitan [13] , Marchenko [3] , and Krein [22] . In geo- 

physics the so-called dynamic deconvolution or layer-stripping (layer-peeling) methods 
emerged, for the identification of layered-earth models from acoustic scattering data 
[2Sl [31j [21 El [5] . More recently the inverse scattering methods have been applied to the 
design and characterization of optical devices involving two interacting modes. Both 
contradirectional coupling and codirectional coupling have been treated. Optical com- 
f— ^ ponents based on contradirectional coupling include thin-film filters and fiber Bragg 

gratings [40l [10l [37l [30l [36l [34] , while codirectional coupling is present in e.g. grating- 
assisted codirectional couplers and long-period gratings [19l [39l [9l [42l [4] . While the 
inverse-scattering problem associated with two interacting modes is well-known, the 
inverse-scattering problem of several, possibly non-dcgcncratc modes (i.e., with dif- 
ferent propagation constants) seems unsolved so far. Some work has been done in the 
case of 4 degenerate modes, that is, two polarization modes in each direction [551 |4"T] . 
and several degenerate modes with only contradirectional coupling pQ. 

On the other hand, several methods for the inverse scattering of acoustic or elec- 
tromagnetic waves in two or three dimensions have been reported. In particular, 
Yagle et al. have developed layer-stripping methods for the multidimensional case 
[45l l43l [44] . By Fourier transforming the problem with respect to the transversal 
coordinates, the multidimensional problem may be regarded as one-dimensional with 
several interacting modes. 
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In this paper we will extend these lines of thought to cover the general inverse scat- 
tering problem associated with any number of interacting modes in one-dimensional, 
reciprocal structures. In the model (Section 2) both codirectional and contradirec- 
tional coupling may be present simultaneously. We limit ourselves to the case where 
the known probing waves and the scattered waves propagate in opposite directions. In 
other words the scattered wave is considered as a reflection from the unknown struc- 
ture. A layer-stripping inverse scattering algorithm is presented in Section 3. Ambi- 
guities related to the simultaneous presence of co- and contradirectional coupling are 
discussed in detail. Possible a priori information that can resolve these ambiguities 
will be suggested. The formalism is particularly useful for the quasi-periodical case 
(Section 4), since only the slowly varying envelope needs to be represented rather 
than the structure itself, yielding an efficient algorithm. In Section 5, the method 
is applied for the numerical reconstruction of a quasi-periodical waveguide structure. 
Sections 4 and 5 are exemplified by a multimode fiber Bragg grating; an optical fiber 
with quasi-periodic refractive index perturbation along the fiber axis, giving rise to 
both co- and contradirectional coupling. Finally, analogies to the multidimensional 
case are discussed in Section 6. 

2. Continuous and discrete coupling model. Consider a structure with P 
modes propagating in each direction along the x-axis. We visualize the x-axis as be- 
ing directed to the right, and say that the +x-direction is the forward direction. The 
propagation constant of the pth. mode is ±/3 p , i.e., the x-dependence of the complex 
field associated with mode p is described by the factor exp(±if3 p x), where the upper 
(lower) sign applies to forward (backward) propagating modes. Note that the prop- 
agation constants of different modes may or may not be different. The propagation 
constants are related to frequency through the dispersion relation of the structure. 
The propagation constants may be expressed p p = n p uj/c, where lo is the angular 
frequency, c is some fixed reference velocity (common for all modes) , and n p accounts 
for the actual phase velocity. (However, in some cases it may rather be convenient to 
express the propagation constants in the form n p uj/c — Tr/A, where A is a constant, see 
SectionHJ) For electromagnetic waves, it is natural to set c equal to the vacuum veloc- 
ity, and consequently we will refer to n p as the effective index associated with mode 
p. In principle, the effective indices may be complex and dependent on frequency, 
meaning that modal loss and dispersion are permitted in the model. However, the 
dispersion must be limited by relativistic causality in the sense that any signal carried 
by the modes travels no faster than the vacuum light velocity. Also, the modal field 
profiles are assumed to have uniform phases such that they can be written real. 

Coupling may occur due to a continuous or discrete scattering structure. In the 
first case, the field is assumed to be governed by the coupled-mode equation 

f=*CE, (2.1) 

where E is a column vector containing the 2P mode amplitudes. In the absence of 
the scattering structure (C CT = C K = 0, see below), the first P elements are the mode 
amplitudes of the forward propagating modes (propagating in the +x direction) and 
the last P elements are those of the backward propagating modes. The coupling 
matrix C can be decomposed into three contributions: 

C = D + C CT + C K . (2.2) 



The contributions can be expressed as 2 x 2 block matrices consisting ofPx? blocks: 



MULTIMODE INVERSE SCATTERING 



3 



D 







(2.3a) 









K* 







(2.3b) 
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(2.3c) 



where * denotes complex conjugate. The first term D describes the frequency depen- 
dence due to the propagation of the different modes ("self-coupling"), and is indepen- 
dent on x; and (3 = diag{/?i, /?2, • ■ ■ , Pp}- Only this term is permitted to be lossy in 
the model (f3 may be complex). In practice, we should require Im/? p L < 1, where L is 
the total length of the structure; otherwise the field at the far end of the structure may 
be close to zero (i.e., the mode will be bound at the left interface to the structure, and 
very little reflection will originate from the far end.). The second term C K describes 
the coupling between counterpropagating modes, whereas the last term C CT accounts 
for the coupling between copropagating modes. The coupling coefficients k and cr are 
dependent on x but assumed independent on frequency. As will become clear shortly, 
the above forms of C K and C CT are consequences of reciprocity and losslessness. It 
should be noted that in structures such as long-period gratings, where the coupling 
is purely codirectional, the coupling is described by n = and a cr with non-zero 
off-diagonal elements. The conventional way of describing such structures would be 
only to consider the upper- left P x P block of C. The layer-stripping method in Sec. 
[3] cannot be used to reconstruct such structures since the reflection response is zero. 

The coupling region in the waveguide is discretized into N layers, each of thickness 
Ax = L/N. If N is sufficiently large so that the matrices in (|2.3p can be treated as 
constants in each layer, we can solve (I2.1j) : 



This transfer matrix relation can be used to propagate the fields through the piecewise 
uniform structure. With the help of the connection between the transfer matrix 
and the scattering matrix (Appendix [B]) we can find the reflection and transmission 
response from the total transfer matrix, obtained as a product of the transfer matrices 
exp (iCAx) of each layer (direct scattering). 

While direct scattering is achieved straightforwardly using the piecewise-uniform 
discretization, for inverse scattering it is convenient to push the discretization further, 
to identify the different contributions to the transfer matrix exp (iCAx). To first order 
in Ax, we have exp(iCAx) = exp(iDAx) exp(iC K Ax) exp(iC CT Aa;). For a continuous 
structure of finite thickness, the bandwidth where the reflection spectrum is signifi- 
cantly different from zero is finite. Thus we need only be concerned with frequencies 
satisfying |o;| < oj\> for some positive constant 0Jb- Note that this model may give 
entirely incorrect results for |a>| > cob- For instance, if P = I the reflection spectrum 
calculated with the discrete model will be periodic with period wc/(niAx), while the 
spectrum associated with a continuous structure tends to zero for large frequencies. 
For inverse scattering, the reflection spectrum and therefore uj^ are known. Therefore, 
provided Ax is chosen sufficiently small we can approximate each layer by a cascade of 
three sections: a section with codirectional coupling, a section with contradirectional 
coupling, and time-delay section. The physical implication of this factorization is that 



E(xj + Ax) = exp (iCAx)E(xj), 



Xj = jAx. 



(2.4) 
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the mode-coupling appears in a discrete point within the layer rather than distributed 
along the whole layer. The contradirectional section may therefore be pictured as a 
discrete reflector. The transfer matrix of the jth layer becomes 



where 

T z = exp(iDAx) = 
T P] = exp(iC K Ax) = 

= exp(iC -Aa;) = 



r 1 o 
o z 

3 1 



*j 
** 



, Z 1 = cxp(i/3Ax) 



(2.5) 



(2.6a) 



*7 



= itanh[(/«*/«) 1 / 2 Aa;](K*/«)- 1 / 2 K* 
tj = cosh[(K*K) 1 / 2 A^]- 1 , 



, $j = exp(ierAx). 



(2.6b) 
(2.6c) 



The form of the matrix in (|2.6bp may for example be verified by evaluating the power 
series expansion of the matrix exponential. In principle, it suffices to express (|2.6[) to 
first order in Ax; however, the exact form is kept to emphasize the properties of each 
of the three sections, to ensure that each section is lossless regardless of the value of 
Ax, and to retain the correspondence to the discrete case (below). 

We are now in the position that we can argue for the forms of the coupling 
matrices (|2.3p . Note that while we have permitted loss in the propagation section 
Z~ x , the coupling sections are assumed lossless. Since the coupling sections also are 
assumed to be reciprocal, their transfer matrices satisfy (|B.I0[) and (jB.ll \ (Appendix 
|B|) . Allowing a more general C re by substituting k* — ► k%i in the (2,1) block, and 
expanding exp(iC K Ax) to first order in Ax, the lossless and reciprocity conditions 
give K12 = — k* and dictate k to be symmetric. Similarly, we can derive the form of 
Co- and establish that «&j must be unitary, i.e., u is hermitian. 

From the discussion above, each layer is characterized by a unitary codirectional 
coupling matrix <f>j and a discrete reflector. Let superscript T denote transpose and let 
|| ■ || be the usual matrix 2-norm. The discrete reflector satisfies Pj = pj and ||pJ| < 1, 
and has an associated, positive definite transmission matrix tj with t 2 = I — PjPj- 

So far we have considered a continuous scattering structure, and discretized it 
into a cascade of codirectional coupling, reflection, and pure propagation. Obviously, 
we can also describe discrete coupling directly. The most general, lossless, reciprocal 
coupling element can be described as a discrete reflector sandwiched between two 
codirectional coupling sections ( Appendix |B|) . Compared to our discrete model above, 
there is an extra codirectional coupling section on the right-hand side of the reflector. 
In the special case where all modes have equal effective index, Z~ x ex J, this coupling 
section commutes with the delay section, and as a result it can be absorbed into the 
next, adjacent layer on the right-hand side. However, in the general case this extra 
coupling section does not commute with the delay section and cannot be ignored. 
For inverse scattering, this coupling section should therefore not be present since 
otherwise, it would not be possible to determine the transmission through the layer 
uniquely from the reflection. Under this assumption, tj is positive semidefinite, and 
uniquely determined by £ 2 = / — pj p* . We restrict ourselves to reflectors that satisfy 
|| p.- 1| < 1; otherwise the reflector will mask the later part of the structure such that 
the inverse scattering procedure will not be possible. Also, with two or more layers 
with \\pj\\ = 1, the structure may behave as an ideal resonator with bound modes. 
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Writing out the transfer matrix (|2.5p of each layer, we obtain 



Z X K, 





ZK 



I 

T, 



— T* 



(2.7) 



where Tj = pj<&j and -RT., = The transfer matrix can be converted into 

a scattering matrix (Appendix [Bj) : 



(2- 



Thus, Tj represents the reflection response from the left of layer j. 

The combined transfer matrix describing the total structure with TV layers is given 

by 



-JV-l J-JV-2 ' ' ' 



Tl Tn 



(2.9) 



From this matrix we can determine the reflection and transmission response using 
(|B.3|) . For example, the reflection response from the left is 



R(u) = S n = -T£t 31 , 



(2.10) 



where Tf-i are the P x P blocks in T. Assuming ||p-|| < 1 for all j, it can be 
proven by induction that T22 is invertible on and above the real frequency axis in the 
complex w-plane, for any number of layers. Physically this is obvious since is the 
transmission response from the right, and therefore it must exist and be causal and 
stable. 

Reciprocity (|B.4a[) gives R(u>) = R(ll>) t . Using ||p-|| < 1 for all j, it can be 
shown by induction that ||it(cd)|| < 1 for a passive structure (a passive structure is 
characterized by lm/3 p > for all p). By causality the reflection response can be 
written in the form 



R(lu) 



h(t) exp(iojt)dt, 



(2.11) 



where h(t) is called the time-domain impulse response. 

When the modes are nondispersive, i.e., (3 is linearly related to frequency, h(t) 
equals a train of non-equally spaced, weighted delta pulses: 



h(t) = ^h k 5{t-t k ). 



(2.12) 



fc=0 



Here h and t k are the weight and arrival time of the kth pulse, respectively. Substi- 
tuting (|2~T2"]) into (|2~TT|) gives 



(2.13) 



fc=0 



The weights h k can in principle be calculated from R(u>) using an inverse transform 
of the form 



h k = 



lim 



1 



♦00 2ui 



max J —uj- 



R{lo) exp(— iujt k )dio. 



(2.14) 
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The arrival times are determined by the delay from a layer to the next of each mode. 
Let At p be the delay of mode p through a single layer. A delta pulse at t = 
is incident to the structure on the left-hand side. Consider the reflection from the 
different layers, as seen from left-hand side of the structure. From layer 0, the arrival 
times in all modes will be zero. An impulse in mode p reflected from layer f into 
mode q, will arrive at At p + At q . Thus, considering layer 1, the arrival times are any 
combinations of two unit delays At p . Considering layer 2, the arrival times are any 
combinations of four unit delays, and so forth. 

When the modes are dispersive, the impulse response is no longer a train of delta 
functions. Nevertheless, for t ~ it can still be written as h°S(t), and the weight h° 
can be found from (|2.14[) . 

Eq. (|2.13p clearly demonstrate that, in principle, for a discrete structure the reflec- 
tion response R(tu) does not approach zero for large frequencies. Only in the special 
case where the modal effective indices are rational numbers with common denomina- 
tor, the reflection spectrum is periodic. Fortunately, in practice, it is not necessary to 
represent the entire bandwidth to enable inverse scattering for a discrete structure. 
As shown in the next section, what is needed in the layer-stripping algorithm is the 
zeroth point of the impulse response, at time t — 0. Since the next nonzero value is for 
t = 2 min p At p Q the zeroth point is computed accurately provided the represented 
bandwidth cj max satisfies w max 3> l/min p At p . Then, if the true reflection spectrum 
is multiplied by a smooth window function W(u>) that goes to zero at u = ±o; m ax, the 
inverse Fourier transform evaluated around zero is approximately w(t)h° , where w(t) 
is the inverse Fourier transform of W(uj). Since 1^(0)^.° « ^ /"™ ax W(aj)R(Lu)duj, 
we can find h° from a measurement of R(lo) in the bandwidth (-w max ,w m ax): 

p»- W{uj)R{uj)dLo 

In many practical cases, the structure to be reconstructed is quasi-sinusoidal. 
More generally, the structure is often quasi-periodic, and e.g. the first "Fourier com- 
ponent" is to be reconstructed. In such cases, one can define modal field envelopes 
which vary slowly with respect to x (compared to a wavelength). Similarly, one can 
extract slowly varying coupling coefficient envelopes. As a result, all quantities in 
(|2.ip vary slowly with x. The relevant bandwidth in (]2 . 1 5[) will then be centered 
about a chosen "design frequency" rather than zero. The main advantage of this 
procedure is that it leads to considerably less requirements on the spatial resolution, 
and as a result efficient inverse scattering. This modification to the model is detailed 
in Section [4] 

3. Layer-stripping method. The inverse scattering problem can now be stated 
as follows: Given a structure consisting of N + 1 layers. Each layer consists of three 
sections (sublayers), the first (3?j) responsible for coupling between copropagating 
modes, the second (pj) responsible for coupling between counterpropagating modes, 
and the third a pure propagating section (Z^ 1 ). The propagation constants of the 
involved modes are known and specified in terms of Z~ . □ From a set of excitation- 
response pairs (that is, R(lo)), we want to reconstruct p.- and &j for all j. 



For simplicity it is assumed a nondispersive structure. 
2 The effective indices may contain small, real, unknown parts An p , i.e., n p = n p i inown + An p , 
where ^ Pi known are known. Provided An p is sufficiently small, the variation of the associated phase 
factor cxp(iLuAn p Ax/c) may be small over the relevant bandwidth. In such cases the unknown parts 
can be absorbed into the #,'s. 
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The structure itself and the medium to the right are assumed to be at rest at time 
t = 0. For incident waves from the left, the reflection response from the structure is 
described by the matrix R(u>) of dimension P x P. This matrix can be viewed as 
the operator which takes the excitation field vector to the reflected field vector. Its 
columns can be interpreted as the responses for orthonormal excitation basis vectors 
ei, e2, . . . , ep, respectively. Here e p has only one nonzero element (equal to unity) 
at position p. Similarly, we can define the forward (uj(uj)) and backward (vj(to)) 
propagating field matrices as P x P matrices where the columns are the fields for 
orthonormal excitations ej., ■ ■ • , ep. A subscript j is specified to emphasize that 
Uj{uj) and Vj(u>) are the fields at the beginning (left-hand side) of layer j. The field 
matrices of layer j + 1 are related to the field matrices of layer j by 



v j+1 (w) 



(3.1) 



where Tj is given by (|2.7p . 

The layer-stripping algorithm is based on the simple fact that the leading edge of 
the impulse response is independent on later parts of the structure due to causality. 
Hence, one can identify the first layer of the structure, and subsequently remove its 
effect using the associated transfer matrix. 

For layer 0, we initialize uo(uj) — I and vq(u>) = R(oj). We define a local reflection 
spectrum Rj(ui) = Vj(ui)uj(uj)~ 1 and the associated impulse response hj(t) as the 
response of the structure after removing the first j—1 layers. Similarly to the impulse 
response of the entire structure, hj(t) contains an isolated delta function at t = 0. 
Due to causality, this pulse is equal to the reflection from the zeroth layer alone. 
Denoting the weight of this pulse Ha, we find from (|2.8[) that 



(3.2) 



Note that Rj{u) is symmetric for all ui as a result of reciprocity; thus h® is symmetric 
as well. Writing out (|3.1[) and (|2.7|) , and substituting Vj{tu) = Rj(uj)uj(uj), we obtain 



u j+ i(u) = Z 1 K j [I - T*Rj(uj)] Uj (u), (3.3a) 

v j+1 (u) = ZK* [Rj(u) - Tj] Uj {uj), (3.3b) 

and therefore 

Rj+i(u) - ZK* [Rj(u) - Tj] [I - T*Rj(lj)] ^ KfZ. (3.4) 

Provided Tj and Kj are known, (|3.4p shows that the local reflection spectrum of layer 
j + 1 can be calculated directly from the local reflection spectrum of layer j without 
calculating the fields itj+i and Vj + \. Note the similarity to the Schur formula used 
in scalar layer-stripping [6]. 

To characterize layer j completely, and to identify Kj , we must determine p ■ and 
By counting the available degrees of freedom (in Tj), we immediately find that 
this cannot be done uniquely. It is therefore necessary to use a priori information on 
Pj and/or &j. The available information may vary from situation to situation. Here 
we will consider the following situations, where Pj and &j can be found using the 
methods in the Appendices I A. 1 1 and IA.21 
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a) = I. In this case there is no codirectional coupling. The identification of 
the layer is now particularly simple, as Pj — Tj uniquely. Note that while 
there is no codirectional coupling, pj describes reflection from all modes into 
all modes. Thus the different modes may still interact. 

b) Pj is diagonal and nonnegative. Now pj is a simple partial reflector which 
only reflects light into the same mode as the incident field (no reflection into 
other modes). The coupling between different modes is instead described 
by <frj. Since Yj = &Jpj&j, pj is found uniquely as the singular value 
matrix associated with Tj, up to reordering of the singular values. Once the 
order of the singular values has been established, the unitary 3>j is found 
uniquely up to the sign of its rows, provided all singular values are distinct 
and nonzero (see Appendix IA.1[) . When one or more singular values of Yj 
are zero, the corresponding row(s) of <&j cannot be determined uniquely. 
More precisely, 3>j is determined up to a premultiplicative unitary matrix J 
operating on the associated mode(s). Physically, this is obvious since when 
a singular value is zero, the associated mode is not reflected from the layer. 
When two or more nonvanishing singular values are equal, &j is determined 
up to a premultiplicative, real unitary J operating on the associated modes. 
Physically, this means that these modes experience the same reflection and 
thus an arbitrary (real) "rotation" of the modes is not detected. In such cases, 
the unitary section &j , as determined by the method in Appendix IA.11 does 
not necessarily correspond to the physical section. This error will propagate 
to the next layers according to (|3.4|) . 

c) &j is symmetric and Pj is real and positive semidefinite. A special case in 
which there are only two degenerate modes in each direction is treated in 
[41]. The reflector matrix Pj can be written PjXjPj, where Pj is a real, 
special unitary matrix and Sj is diagonal and nonnegative. Since Tj = 
<&Jpj<&j = $JPjSjPj$j, we find £j and Pj<&j as in the previous case, 
with the identical ambiguity issues. The separate identification of Pj and <&j 
is accomplished using the factorization method in Appendix IA.2( with certain 
ambiguities related to the sign of the eigenvalues of 

The ambiguities when determining &j in situation b) are in fact very similar to 
the well-known ambiguities in the scalar case with a single mode in each direction. In 
the scalar case any ir phase-shift sections between the reflectors cannot be identified 
since the associated round-trip phase accumulated to and from a reflector becomes 
2-7T. In our multimode case, the sign of the rows of the "phase-delay" section (4>j) 
between two reflectors cannot be identified. Similarly, in the scalar case, any phase- 
shift section preceeding a zero reflector cannot be determined uniquely. Instead it is 
chosen arbitrarily (e.g. removed), and attributed to the next layer with a nonzero 
reflector. 

When the structure to be reconstructed is a discretized version of a smooth struc- 
ture, the smoothness can be used to resolve ambiguites. First we consider situation 
b). For small Ax, <frj is close to identity; thus the sign of the rows of 3>j can be 
determined uniquely. If Pj has distinct eigenvalues, valid for all j, the order of the 
eigenvalues of p can be determined from the order of the eigenvalues of Pj_\ using 
the smoothness of k = k(x). If there are equal eigenvalues for a certain reflector p., 
or if Pj is singular, the ambiguities of *f>j are characterized by the premultiplicative J 
matrix (Appendix lA.l[) . In other words, the chosen <f>j is related to the corresponding 
true matrix (^f^true) by 4?-, = J<&j.tmc- By choosing J such that — 3>j-i|| is min- 
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imum, the resulting J is close to identity (that is, || J — 1\\ < 2||<&j true — <&,-_x||). Since 
tj and Z -1 are close to identity as well, the order of three sections «7, tj, and Z~ 
can be interchanged (see Section [5]) . Thus the error due to wrong choice of &j can 
be absorbed into ■ More generally, provided only a few neighboring layers have 
singular or degenerate Pj's, only the corresponding and following &j sections may be 
determined erroneously, and the determination of the later part of the structure is 
(approximately) unaffected. 

In situation c), the order of eigenvalues of pj can be determined as in situation b). 
However, Pj&j is not necessarily close to identity. Nevertheless, the sign of its rows 
can be determined from if k = k(x) is sufficiently smooth. (Recall that 

Pj&j is unitary, which means that in each row there exists at least one element of 
magnitude > 1/yT.) Finally, since 3>j is close to identity, its eigenvalues are close to 
unity. It follows that the factorization of Pj&j into Pj and 3>j is unique (Appendix 

E3- 

From the discussion above, we summarize the layer-stripping algorithm, analo- 
gously to the scalar version described in ref. [6l [5] , that can be applied to identify a 
structure supporting multiple modes: 

1) Initialize j = 0. Set Rj{u>) — R(lu). 

2) Compute the zeroth weight h® of the impulse response. In practice this is 
achieved by the substitutions h° — > h® and R(uj) — > Rj(u) in (|2.15[) . 

3) Use a model-specific factorization of h° = 3>J Pj&j to find <f?j and p -. 

4) Calculate tj = (I — PjP*) 1 / 2 such that the associated eigenvalues are positive, 
and set Kj = tj <&j. 

5) Calculate the next, local reflection response Rj + i(uj) using (|3.4p . 

6) If j < N — 1, increase j and return to [2] 

When the scattering structure is continuous, one can use the true reflection spec- 
trum as input to the layer-stripping algorithm, even though the structure is modelled 
discrete. This can be justified as follows: The layer thickness Ax is chosen small such 
that the first order approximations of exp(iC K Aa;) and exp(iC CT Ax) are accurate. 
(Thus an upper bound on ||C K |j and ||C CT || should be known a priori.) Let u < lo^ 
be the bandwidth where the true reflection spectrum is significantly different from 
zero. For sufficiently small Ax, the first order approximation of exp(iDAx) is valid, 
and the true reflection spectrum is approximately equal to that of the corresponding 
discrete model in the bandwidth w < Wb- In the limit t — > + , the (p, q) element of 
the impulse response of the continuous structure can be calculated exactly from (|2.1[) 
using the Born approximation, yielding 

1 f°° 

h pq {t = + ) = — lim / R pq (uj) exp(— iuit)duj = in* (x = + )c/(n p + n„). (3.5) 

27T (-.0+ J^oo 

Here K pq (x — + ) is the (p, q) element of k(x) at x = + . For practical computations, 
the integral in Eq. (|3.5[1 must be truncated at ±Wb! thus, to find the leading edge of 
h pq (t), one can take t = in the integral, and multiply the result by a factor of two. 
(Recall that by causality lim t ^ - h pq (t) = 0.) Once k for the zeroth layer is found, 
one can propagate the fields using (|3.4p . Since we have not identified the codirectional 
coupling $o of the zeroth layer, $o is associated with the next layer. Thus, after the 
zeroth layer has been stripped off, the leading edge of the impulse response of the 
remaining structure becomes 

*J [Kq( x = Az + )c/(n p + n q )\ *„, (3.6) 
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where the square bracket denotes a matrix formed by the elements inside. The identi- 
fication of 4> an d [iKp q /(n p + n q )~\ can now be accomplished using the factorization 
methods described above. The algorithm continues in the same way, until finally the 
bandwidth of the reflection spectrum of the remaining structure exceeds This 
remaining part of the structure can be made arbitrarily thin by choosing a sufficiently 
small Ax. 

The difference between the latter "quasi-continuous" formulation and the discrete 
algorithm is essentially the factor n p + n q , and the method for evaluating the leading 
edge or first point of the impulse response. When the effective indices can be approxi- 
mated by some number no for allp, n p ~ no, one can in fact use the discrete algorithm 
directly: A periodic extension of the true reflection spectrum outside a principal band- 
width [— Wmax.Wmax] corresponds then to a discrete model with Ax = 7rc/(2now m ax). 
The first point of the impulse response is calculated by (|2.15j) using a rectangular 
window function W(lo). For a broad class of waveguides of practical interest, the 
effective indices are similar (see Section 4). While the phase relation between the 
modes, as described by Z~ , may still result in a nontrivial multimode coupling, the 
discrete algorithm gives accurate results. The errors due to this periodic spectrum 
approximation can be corrected to some extent by including the factor (n p +n q ) / (2n ) 
in the elements on the right-hand side of (|2.15|) . This can be justified e.g. using the 
Born approximation. 

4. Quasi-sinusoidal coupling structures. Continuous coupling in acoustical, 
radio frequency, or optical waveguides may be obtained by perturbation of the effective 
indices n p associated with each mode. This can be achieved by modulation of the 
wall profile or waveguide medium properties. As a concrete example, we will discuss 
fiber Bragg gratings [17] . which have attracted large interest recently due to their 
applications in fiber optical communications and sensors. A fiber grating is formed in 
an optical fiber by modulating the refractive index of the core periodically or quasi- 
periodically. The main peak of the reflection spectrum appears for the frequency where 
the reflection from a crest in the index modulation is in phase with the next reflection. 
Permanent gratings are fabricated by UV-illumination. In fibers doped with certain 
dopants such as germanium, the UV-illumination will permanently rise the refractive 
index of the core. Advanced fabrication methods have made it possible to manufacture 
complex gratings with varying index modulation amplitude and period. The layer- 
stripping algorithm is the most widely used method for designing the index profile to 
obtain a given reflection spectrum [lOj [37l [36] . 

In most cases, the fiber grating is formed in a single-mode fiber, and coupling is 
only considered between the forward-propagating and backward-propagating funda- 
mental mode. The field matrices Uj(u>) and Vj(u>) are then scalar functions. However, 
in some cases it is not sufficient to consider only one forward-propagating mode and 
one backward-propagating mode. For instance, a single mode fiber is always slightly 
bircfringent, and the photosensitivity can be polarization-dependent [16 . In this case, 
two forward-propagating and two backward-propagating polarization modes must be 
considered. An inverse scattering algorithm that takes into account polarization mode 
coupling is described in |41j . The coupling between the two polarization modes are 
described by Jones matrices |20) . Both polarization modes have approximately the 
same effective index, so Z~ x — exp(i(3Ax)I , where the common propagation constant 
P is scalar. 

In a multi-mode fiber, the modulation of the refractive index may result in cou- 
pling between the fundamental mode and other modes. Each mode has a transversal 
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field profile ^ p {r,6) which is a solution to the scalar wave equation in polar coordi- 
nates r and 4> [38j3 

{ V? + k 2 (n 2 (r) - n 2 p )} 9 p {r, <f>) = 0. (4.1) 

Here n(r) is the unperturbed, refractive index profile of the fiber, which is assumed to 
be real, Vt is the transversal nabla operator, and k — lu/c. The field ^ P {r, <f>) and its 
first derivatives are continuous. For bound modes, the fields are real and orthonormal 
such that J A *& p (r, cf>)^ q (r, <f>)dA — S(p — q), where 5(p — q) denotes the Kronecker 
delta, and Aoa is the entire transversal plane. The effective indices n p are eigenvalue 
solutions to (|4.ip . A mode p is bound when n c \ < n p < n co , where n co and n c \ are 
the refractive indices of the fiber core and cladding, respectively. Ignoring radiation 
modes, which in the vicinity of the core decay rapidly away from the excitation source, 
the total electric field E(r, <fi, x) can be written as a superposition of forward- and 
backward-propagating bound modes: 

p 

E(r,cf>,x) = 5>+(ar) + b p (x))* p (r,cf>). (4.2) 
P =i 

Here b p (x) contain all x-dependence including the harmonic propagation factor exp(± i/3 p 
where f3 p = kn p . 

Coupling between the modes originates from longitudinal modulation of the re- 
fractive index. Let the refractive index be perturbed quasi-periodically with a spatial 
period A, 

n(r, 4>, x) = n(r) + An ac (r, </>, x) cos + S(x)^j + An dc (r, (f>, x), (4.3) 

where An ac (r, (f>, x), An<i c (r, <f>, x), and 8{x) are slowly varying with x over a distance 
A. We assume that An. ac (r, <j>, x) n, and Artdc(^, 0, x) -C n, which is the case for 
practical fiber gratings. The total electric field must satisfy the scalar wave-equation 
for the perturbed fiber, i.e., 

V ? + ^2 + fc V ( r > <t>> x ) } E (r, 4, x) = 0. (4.4) 

We now substitute (|4. 2|) into (]4.4[) . take (|4. 1 [) into account, and multiply the resulting 
equation by ^ q (r, <j)). By integration over the entire transversal plane, and recalling 
that the modes are orthonormal, the resulting set of second order differential equations 
can be decomposed into first order coupled mode equations |38j . 



M+(x) 
dx 



iP P b+(x) =iJ2c pq (x)(b+(x) + b q (x)), (4.5a) 



9=1 
P 



^^)- + p b-(x) = -iJ2C pq (x)(b+(x)+b-(x)), (4.5b) 



9=1 



3 To find the exact electromagnetic modes, the vector wave equation must be solved. However, 
for weakly guiding waveguides (waveguides with small difference between the refractive index of the 
core and the cladding), the scalar wave equation can be used. This is the case for most conventional 
fibers. 



12 



O. H. WAAGAARD AND J. SKAAR 



where 

C Pq {x) = f (n 2 (r, <t>, x) - n 2 (r))^ p (r, cf>)* q (r, </>)dA. (4.6) 

Note that the frequency-dependence of (|4.6p can be ignored in practice, since the 
normalized bandwidth of interest is usually much less than unity, and the field profiles 
and effective indices are approximately constant in this bandwidth. Also note that 
since the fiber is assumed to be weakly guiding, n p can be set equal to n co ; thus 
C — C 

In the case of a quasi-periodic structure it is natural to write the coupling coeffi- 
cient as a quasi-Fourier series: 

2vr \ _ , / 2tt 



Cpq(x) = (Jpq{x) + Kp q (x) exp ( i~^X J + K* q (x) 6Xp I -i~£ X 



J2 ^^exp 

m|>2 



lirm 
i — ; — x 
A 



(4.7) 



where the "Fourier coefficients" K pq (x), a pq (x), and Kpq\x) are slowly varying over 
a period A. For a fiber grating the index modulation n(r, <f>, x) — n(r) is given by 
(|4.3|) and is small compared to n(r), so the zeroth and first order Fourier components 
dominate. Note that &rg{K pq (x)} — 9(x). 

The field amplitudes b p (x) vary rapidly; it is therefore convenient to introduce 
the slowly varying field envelopes u p {x) and v p {x) by setting 



<+{x) = i 1/2 u p (x) exp (ijx^j exp , (4.8a) 

r(x) = r 1/2 v p (x) exp (-» Js) exp f-i^T\ . (4.8b) 



Since an identical phase factor is removed from all modes, the reflection response as 
calculated from b p and b~ will only differ from that calculated from u p and v q by 
a constant phase factor not dependent on p and q. Inserting (14. 7|) and (|4.8p into 
(|4.5p . and ignoring rapidly oscillating terms (since they contribute little to du p /dx 
and dv p /dx), we obtain an alternative set of coupled- mode equations 

= i5pU p {x) - ^~~^~ u p( x ) +i^2^ P q{x)u q (x) +y^\K pq (x)\v q (x), (4.9a) 

q=l q=l 

= -id p v p (x) + \—^ v p( x ) -i^Vpqi^Vqix) +^2\ K pq(x)\u q {x), (4.9b) 

q=l q=l 

where S p — (3 P — tt/A = n p uj/c — n/ A is the wavenumber detuning of mode p. Thus, 
— ?|fCp g (iE)| is the coupling coefficient between modes p and q propagating in opposite 
directions, while <J pq (x) — 8{p — q)(d0(x)/dx)/2 is the coupling coefficient between 
modes p and q in the same direction. With E = [ui, U2, . . . ,Up, U2, • • • , vp] T we 
find that ()4.9|) coincides with (|2.ip . where <r pq {x) — 8{p— q)(d0(x)/dx)/2 and —i\K pq (x)\ 
are the (p, q) elements of cr and k, respectively, and S p are the diagonal elements of 
(3. Note that S p do not correspond to the actual propagation constants but rather 
their detuning from w/ A. Approximating the effective indices by n co , this means that 
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the bandwidth of interest is not centered about zero but rather about the "design 
frequency" oj = nc/(n co A). The frequency interval of integration in (|2.15p should 
be centered about ujq. As in the scalar case we also note that in general, the 
geometrical phase variation 6{x) cannot be distinguished from the phase variation 
associated with the dc index term An<jc(?", <fi, %)■ 

We observe that er is real and symmetric, and k is imaginary and symmetric. 
Moreover, it is not difficult to realize that ik is positive semidefinitecl Thus <J>j defined 
in (|2.6p is unitary and symmetric, and —pj is real and positive semidefinite. It follows 
that we can use the layer-stripping method together with the factorization approach 
c), as given in Section^ to identify the coupling sections p^ and 3>j (and therefore the 
coupling matrices k and er as a function of position x). Since (i<&j) T (— Pj)(i&j) = 
<&J Pj<&j, the factorization approach gives — Pj and i&j. 

For a fiber grating it is usually reasonable to assume that the ac and dc in- 
dex modulations can be written in the forms An ac (r, 0, x) = An(r, 4>)An ac (x) and 
Andc(^, 0, x) = An(r, <p)Andc(x), respectively. Here An(r, (f>) accounts for the transver- 
sal variation of the index modulation profile, and An ac (a;) and And c {x) are the ac and 
dc modulations as a function of x. As before, we assume that the index modulation 
and n co — n c \ are small, yielding 

.An ac (x) , . 

= -* 2 Vl ^ ^ 

ct(x) = An dc (x) V - (4.10b) 
2 dec 

where rj is independent on x. The elements of rj are 

V Pq = k An(r,(f))%^ q dA. (4.11) 

When the mode profiles and An{r, <p) are known, this means that the entire coupling 
matrix k(x) is determined from only a single nonvanishing element. For er, two 
elements are needed (including at least one diagonal element). Note that in this case, 
it is indeed possible to distinguish between the dc index modulation An<i c (x) and the 
geometrical phase variation d9(x)/dx using information contained in er. 

For characterization of multimode gratings, measurements of the reflection from 
every mode to every mode are required. Performing such measurements is not triv- 
ial. In Ref. |33j . an auxiliary long-period grating (LPG), i.e, a grating with purely 
codirectional coupling, is used to characterize another interrogated LPG. Fig. 14.11 
shows how this method can be adopted to characterization of multimode fiber Bragg 
gratings (FBGs) using optical frequency domain refiectometry 12J, provided there are 
no degenerated modes. Light is coupled into the fundamental mode of the input fiber 
and the frequency of the highly coherent source is swept. The coupler splits the light 
equally into two fibers. The LPG couples light from the fundamental modes into the 
other modes so that the total optical power is distributed between all modes. The 
light returned by the FBG will again propagate through the LPG, and some light 
from each mode will be coupled back into the fundamental mode. The mirror reflects 
only the fundamental mode, and at the coupler the reflected light from the mirror 

4 The real matrix given by the elements ^p^q is clearly positive semidefinite, since 
Sp, q a p*p*9 a 9 = (Ep *p"p) 2 > for any real a p . For a fiber grating An ac (r,ijl,3;) > for 
all r and <f>\ thus \K pq (x)\ adopts the positive semidefinite property from fyptyq. 
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Coupler 




Fig. 4.1 



Measurement setup for characterization of multimode gratings. 



interferes with the light in the fundamental mode out of the LPG. If the fiber between 
the LPG and FBG is sufficiently long such that the difference in delay between the 
modes is larger than the length of the impulse response of the FBG, the individual 
elements of the reflection matrix will be separable in the time-domain. 

5. Numerical example. A potential application of the multimode layer-stripping 
method is to characterize coupling from the core mode to cladding modes in a single 
mode fiber. Cladding modes are not bound within the core of the fiber, but by the 
cladding/air boundary [8]. A single mode fiber may support as many as 100 cladding 
modes. The power in these modes will eventually be lost to the environment. The 
core-cladding mode coupling can be seen clearly in the transmission spectra of strong 
gratings. For chirped gratings [26] and chirped, sampled gratings [27], the bandwidth 
may become larger than the separation in resonant wavelength between the core-core 
mode coupling and the core-cladding mode coupling. Then the core-cladding mode 
coupling will interfere with the reflection spectrum associated with the core mode 
[llj . This unwanted coupling is often handled by writing the grating in fibers with 
depressed cladding modes [7]. There has also been some attempts of taking into ac- 
count the core-cladding mode coupling in the design of the grating [23] [14]. Here, 
direct scattering is treated with multiple mode coupling, but the inverse scattering 
has so far been purely single-mode. The layer-stripping algorithm described in Sec- 
tion [3] can be used for characterization of such coupling and possibly for design. In 
contrast to the methods in [23] [14] , multiple modes can be taken into account in the 
inverse scattering part of an iterative design process. 

A simpler, but nevertheless interesting problem is to characterize coupling in an 
optical fiber with a few bound modes. Here, we will present a numerical experiment 
simulating a grating in a fiber with n co = 1.452, n c \ — 1.437, and core radius r co = 
5 /im. By solving the eigenvalue equation for a circular fiber [38] , we find that this 
fiber supports four modes: LP i, LPn, LP 2 i, and LP 2 at the design wavelength Ao = 
1.55 /im. Here, the index I in LP/ m means that the transversal field profile can be 
written in the form /z m (r) cos(7</>). In the further discussion, these modes are denoted 
1 to 4 in the order indicated above. The eigenvalue equation gives the modal indices 
rti = 1.449, 7i2 = 1.444, ri3=1.439, and ri4=1.437. We assume that the refractive index 
is modulated uniformly in the core of the fiber, but not at all in the cladding. This 
is quite realistic since, during fabrication, the fiber usually is made sensitive to UV 
exposure only in the core. By evaluating (|4.11[) , we find that there will be no coupling 
between modes with different azimuthal indices I: 

0.957 -0.116" 

0.874 , , 

0.707 ' ^ ' 

-0.116 0.491 

There is no coupling to or from modes 2 and 3; thus the grating profile can be found by 
applying a scalar layer-stripping method separately to the responses associated with 



2tt 
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these modes. On the other hand, modes 1 and 4 are coupled, so that the multimode 
layer-stripping method must be applied when using the associated responses as a 
starting point. 

Defining the nominal mode index uq = (ni + n 4 )/2, the grating period is set to 
A = Ao/(2no). The length of the grating is L = 20 mm, and An ac (a;) has the form of 
a raised cosine window with maximum value 1 • 10~ 3 . Furthermore, An<ic(x) is chosen 
as a sine-modulated Gaussian window with full-width-at-half-maximum of 7 mm and 
a maximum value 5 • 10 -4 ; the period of the sine-modulation is 4 mm. The grating is 
chirped by varying the grating phase according to 



The reflection matrix as a function of frequency detuning is generated using the piece- 
wise uniform approximation (Section^) with Ax = 10/^m, which gives ./V=2000. Zero 
detuning is taken to be the frequency /o = c/\ . Figure IBTTk ) shows the resulting re- 
flection matrix spectrum. The maximum values are [|i?n|, I-R22I, | -R33 ] , I-R44I, |-Ri4|]max 
[99.6, 99.6, 97.0, 83.0, 28.3]%. Note that the large chirp has resulted in significant 
spectral overlap between the different elements. 

The reflection matrix is applied as input to the layer-stripping method. As the 
modal indices are similar in magnitude, we use the discrete algorithm directly, and 
Yj is calculated by taking into account the factor (n p + n q )/(2n ) as discussed in 
Section [3J Moreover, n(x) and <r(x) is calculated by inverting the expressions for pj 
and 3?j in (|2.6b[) and (|2.6cp . respectively. Figure l5~Tb ) shows An ac (x) along with its 
reconstructed version. The reconstructed An ac (a;) is calculated by a least square fit to 
(|4.10ap using the diagonal elements of the reconstructed k(x). We find that the error 
in reconstructed profile is less that 4 • 10~ 6 m _1 . Also shown is the ac modulation 
profile calculated using scalar layer-stripping on Rn. Due to the strong coupling 
between mode 1 and 4, the scalar layer-stripping method does not reconstruct the 
profile accurately. Figure 15.1b ) and I5.1b l) show that it is possible to separate the dc 
index variations Andc(x) from the grating phase gradient d6(x)/dx. The separation 
is based on a least square fit to (|4.10bj) using the diagonal elements of er{x). The error 
in reconstructed An& c (x) is less than 6 • 10~ 5 m , while the error in reconstructed 
A9(x)/Ax is less than 300 m _1 . Errors are mainly due to the finite Ax in addition 
to the fact that the reflection matrix spectrum of the discretized structure is strictly 
nonperiodic (see last paragraph of Section [3]) . 

6. Analogies to 3D inverse scattering. An important inverse scattering 
problem is the three-dimensional problem associated with the Schrodinger equation 



where ip(x,y,z,k) is the wave function and V(x,y, z) is a smooth and nonnegative 
potential with compact support. In particular, solutions to this problem is applica- 
ble to inverse seismic scattering. This problem has been solved using a generalized 
Marchenko method in [55] and [32], while layer-stripping solutions are suggested in 
[4"5] and [33]. Note the close resemblance between (|6.1I) and (14.41) . indicating that a 
similar method as that in Section [4] can be used. 

We express the solution as a superposition of the eigenmodes of the Schrodinger 
equation with V(x,y, z) — 0. Writing ip(x,y, z;k) = ^(y, z; k y , k z ) exp(ik x x), these 




(5.2) 



{V 2 + k 2 -V(x,y,z)}^{x, y> z;k) = 0, 



(6.1) 
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Fig. 5.1. a) Magnitude of the reflection spectrum \Rn\ (solid curve), |i?22| (dashed curve), 
I-R33I (dashed- dotted curve), | -R44 1 (solid point-marked curve), \Ri4\ = \Ra\ (dotted curve), b) Re- 
constructed longitudinal ac modulation An a c(x) (solid curve), actual ac modulation (solid point- 
marked curve) and ac modulation calculated using scalar layer- stripping on Rn (dashed- dotted 
curve), c) Reconstructed longitudinal dc modulation An,j c (i) (solid curve), actual dc modulation 
(solid point-marked curve) d) Reconstructed grating phase gradient d6 / dx (solid curve) and actual 
grating phase gradient (solid point-marked curve). 



eigenmodes are given by 

*(y, z; k y , k z ) = exp(i(k y y + k z z)), (6.2) 

where k y and k z are the wave numbers in y-direction and z-direction, respectively, 
and k 2 = k 2 x + kl + k 2 . 

In a discrete model, the wavenumbers k y and k z can for example be discretized in 
equal intervals Ak, such that k y = pAk and k z = qAk. In the j/z-plane, this means 
that only a principal range (— ir/Ak, -k / Ak) is considered, and the fields are extended 
periodically outside this range. The integers p and q are the modal indices satisfying 
p 2 + q 2 < (k/ Ak) 2 for propagating (not evanescent) modes. The modal field profiles 
are written in normalized form ^ pq {y 1 z) = (Afc/27r)^I'(j/, z;pAk, qAk). The total field 
ip(x, y, z; k) is expressed as the superposition 

i>(x, y, z; k) = J2( b pg ( x ) + b p^ x ))^Pi^^ z )> ( 6 - 3 ) 
p-.q 
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where b p t q (x) includes all x-dependence of the fields, and ± indicate the sign of k x , 
i.e, the propagation direction of the mode. 

As in Section HI we insert (16. 3p into (|6.1|) . multiply by ^* g (y,z) and integrate 
over the principal range of the yz-plane. This leads to the coupled mode equations 

i^C pq>ra (x)(b+(x) + b-{x)), (6.4a) 

r,s 

(x)+b-Jx)), (6.4b) 

r,s 



dx 



^x,pqbpg{x^j — 



^ x ' <' n >x,pq u pq\ u 'J 



where the coupling coefficients are given by 

C pq , rs (x) = -^r / % q (y,z)V(x,y,z)^ rs (y,z)dydz 



1 / Ak\ 2 f 

~2k~ \2tt) J V ( x ' y > z ^ ex v{ iAk (( r ~ p ' >y+ ( s ~ q ') z ^ dydz > 



(6.5) 



and fc r 



= [k 2 — (Ak) 2 (p 2 + q 2 )] 1 / 2 . We restrict ourselves to the situation where 
V(x,y,z) is known to be quasi-periodic along the x-direction. Then an expansion 
of C pq<rs {x) as in (|4.7j) together with the transformation (|4.8p can be used, resulting 
in the exact same problem as that described in Section 2] Thus the layer-stripping 
method in Section [3] can be applied. The required input data is the reflection into all 
plane waves upon excitation of the different plane waves onto the plane x = 0. The 
scattering potential V(x, y, z) is found from the inverse of (|6.5[) . 

There are two complications. First, in order to use the factorization methods 
developed in Section [3l we must ensure that reciprocity implies symmetric scattering 
matrices. This is guaranteed when the mode profiles can be written real. Thus we 
define real mode fields by the transformation 



(6.6) 
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++ denotes a column vector containing the modal field amplitudes ^ pq with 



Here 
positive p and q 



\1/ |_ denotes a column vector containing the modal field amplitudes 



with negative p and positive q, and so forth. The dimension of the identity matrices in 
the blocks of M corresponds to the dimension of If C denotes the matrix formed 

by the elements C pq , rs , the coupling matrix transforms C — > M*CM T . Inspection of 
(|6.5p shows that the transformed — C is real and positive semidefinite (recall that 
V(x,y,z) > 0); thus enabling the factorization method in Section[3l 

Second, the causality argument of the layer-stripping method does only work 
when the coupling matrix C is independent on frequency. Eq. (|6.5|) shows that this 
condition can only be justified when the relevant frequency band is narrow. Therefore 
the structure must, in addition to be quasi-periodic along the x-direction, vary slowly 
along the transversal direction. The variation must be sufficiently slow such that the 
modes with (p 2 + q 2 )Ak 2 <C k 2 contain sufficient information about the transversal 
dependence, and the other modes may be neglected. 
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7. Conclusion. A layer-stripping method for the inverse scattering of multi- 
mode structures has been proposed. Ambiguities related to factorization of each 
layer's response into codirectional and contradirectional coupling have been discussed. 
When there are no codirectional coupling, the ambiguities disappear. Also, when the 
structure to be reconstructed is smooth, there are important cases with simultaneous 
co- and contradirectional coupling that can be reconstructed uniquely, provided the 
reflector eigenvalues are nonzero and nondegenerate. Applications to quasi-periodical 
structures, and analogies to multidimensional inverse scattering have been discussed. 

Appendix A. Matrix factorizations. 

A.l. Takagi factorization of complex symmetric matrices. Any complex 
symmtric matrix T can be written 



Y = U L YXJ, 



(A.l) 



where U is unitary and X is diagonal and nonnegative (See e.g. |18) . Chapter 4.4). 
Eq. (jA.ip is called Takagi factorization. 

A constructive proof, suitable for implementation, can be given as follows: Sin- 
gular value decomposition yields 



T = V{E,V 2 



(A.2) 



Y T and 



where Vi 2 are unitary, and S is diagonal and nonnegative. Using Y 
(TY t ) T = Y f Y we find that WT, = T,W T = SW, where W = V* 2 V x . Thus, pro- 
vided Y is nonsingular, W is symmetric. Then VW can be chosen such that it com- 



mutes with X and is symmetric, and we obtain Y = 
or 

Y = U T HU, 



V^WEV 2 = (VWV 2 ) l 2~:VWV 2 , 



(A.3) 



where U = \/WV 2 is unitary and S is diagonal and positive. 
If Y is singular, we write 







and 



W 



Wn w 12 
W 2 i w 22 



(A.4) 



where we have arranged X so that the zero singular values are the last ones, S' is a 
diagonal matrix with the nonzero singular values, and Wn has the same dimension 
as We now find E'W n = W n £', W 12 = W 21 = 0, and W n = Wj v The 
commutation relations do not provide any information on W 22 - Choose \JW such 
that 







Vw^ 



(A.5) 



where y/Wu is symmetric and \fW\\ and S' commute. Write Y = UlY,U 2 , with 

U 1 = VW T V 2 = ^ (A.6) 
U 2 = VWV 2 = \ ¥ T „ ] . (A.7) 
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The matrices U'[ and U 2 are the rows of TJ \ and U 2 that correspond to the zero 
singular values, and they do not give any contribution to T. We may therefore 
replace the rows U" by U'2, which gives U\ = U 2 = U. 

The matrix S is unique up to reordering of the singular values. When the order of 
the singular values is established, U is unique up to the replacement JU — > U, where 
J is a unitary matrix satisfying (JU) T HJU = C/ T X£7. This leads to J T SJ = E. 
Assuming the singular values are sorted in, say, descending order, we find that J is a 
unitary block-diagonal matrix, where each block has a dimension equal to the number 
of corresponding repeated singular values. For zero singular values, the corresponding 
block in J is an arbitrary unitary matrix. For repeated non-zero singular values, 
the corresponding block in J is real. For a distinct, non-zero singular value, the 
corresponding block of J is either 1 or — 1. 

A. 2. Factorization of a unitary matrix into a symmetric matrix and an 
orthogonal matrix. A unitary matrix U can be factorized into U = where 
P is a real unitary matrix (orthogonal matrix) and $ is a symmetric unitary matrix 
(See e.g. [18] . Chapter 3.4). A constructive proof, suitable for implementation, can be 
given as follows. First we note that the symmetric unitary matrix $ can be factorized 
into = P\DP\, where D is a diagonal unitary matrix and P± is a real unitary 
matrix (a simple, constructive proof for this particular spectral decomposition is given 
in [TS], Chapter 4.4). Thus, an equivalent problem is to show that 

U = P 2 DPj, (A.8) 

where P2 = PPi- The decomposition in (|A.8[) is very similar to singular value 
decomposition of real matrices, except that D may have complex elements. 
The matrix U T U is unitary and symmetric; thus we can write 

U T U = PiAPj, (A.9) 

where P\ is a real unitary matrix and A is a diagonal unitary matrix. Define 

P 2 = UP 1 D*, (A.IO) 

where the diagonal matrix D is a solution to D 2 = A. The matrix P 2 is unitary since 
it is produced by multiplication of unitary matrices, thus P^Pj = I- The matrix is 
also real since 

P2P2 = D*PjU T UPiD* = D*P^P 1 D 2 PjP 1 D* = I, (All) 

which giv es P 2 = (P;Pj)P 2 = P;(PjP 2 ) = P* 2 - 

From (jA.lOp we therefore conclude that the decomposition (|A.8|1 . with real unitary 
Pi and P2 and diagonal _D, is always possible. It follows that any unitary matrix can 
be written U = P<I>, where P is real and unitary, and $ is symmetric and unitary. 
Note that any global phase of P can instead be assigned to so without loss of 
generality we can assume that P is special (det P = 1 and det $ = det U) . 

Since D is calculated from D 2 = A, the sign of its elements are arbitrary. The 
ambiguities when determining P\ in (|A.9[) give rise to ambiguities in P and The 
possible P and $ can be expressed as P = UP±JD* J T Pj and $ = P\3DJ T P\ 
for a real unitary J that commutes with D 2 . Here P\ is fixed. If the signs of the 
elements of D are known to be such that any equal elements of D 2 correspond to 
equal elements of _D, then J commutes with D and can be ignored. 

Appendix B. Linear, reciprocal and lossless components. 
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Consider a linear component with P input and P output modes on the left-hand 
side, and also P input and P output modes on the right-hand side, see Fig. IB. II The 



s 



Fig. B.l. A linear component with P input and P output modes on each side. 
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component is completely characterized by the IP x 2P dimensional scattering matrix 
S which relates the input and output fields: 



(B.l) 



The field vectors that propagate to the right and left are denoted u and v, respectively, 
and the subscripts 1 and 2 indicate the left- and right-hand side of the component. 
The scattering matrix is a block matrix; the blocks Sn and S22 being the reflection 
from the left and right side of the device, respectively, and S21 and S12 the transmis- 
sion through the device from the left and right, respectively. These blocks have the 
dimension P x P. 

There exists a similar relation, a transfer matrix relation, that connects the fields 
on the left-hand side to the fields on the right-hand side: 



U2 


= T 


Ui 


v 2 _ 




_Vl_ 



(B.2) 



(B.3) 



Comparing (|B.1[) and (|B.2j) wc find the blocks of T: 

£21 — £'22>S : i2 L 'S'll S22S12 1 

— S 12 Sn S 12 

To describe a device with a transfer matrix, S12 must be invertible, that is the trans- 
mission from the right cannot be zero for any input field vector. Thus, ideal mirrors, 
for example, cannot be described by a transfer matrix. 

Provided the mode profiles can be written real, reciprocity means that the scat- 
tering matrix is symmetric [29l 115) . i.e., 

£11 = 



(B.4a) 



S22 — S 



(B.4b) 

5 2 i = Sj 2 . (B.4c) 
Moreover, the lossless condition is expressed as the unitarity condition S^S = I: 

Sl 1 S u +Sl 1 S 2 i=I (B.5a) 
S\ 2 S 12 + Sl 2 S 22 = I (B.5b) 
S\ 2 S n + Sl 2 S2i=0. (B.5c) 
With (|B.4[) in mind, we introduce Takagi factorization of Sn and —S22 (see Appendix 



S22 = M-p')®! 

S21 = $ r t'*i. 



(B.6a) 

(B.6b) 
(B.6c) 
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Here, 4>i and <fr r are unitary matrices, p and p' are diagonal and nonnegative, and 
t' = <&JS 2 i*i- By substituting into (|B.5|) and using (|B.4|) we obtain 

t'V = 1 p 2 (B.7a) 
t't't = 1 p' 2 (B.7b) 
p' = t'pt'*- 1 . (B.7c) 

Introducing the singular value decomposition i' = C/'iV', we obtain from (|B.7a[) that 
t 2 = V'(I - p 2 )V'\ which means t = V' y/I - p 2 V n . Backsubstitution shows that 
t' can be written t' = V 1 — p 2 for a unitary U; thus (|B.7c[) reduces to p' = UpU T . 
With these properties, it is straightforward to show that (|B.6|) can be written 

S u = <f>Jp$i (B.8a) 

S 2 2 = *r(-p)*^ (B.8b) 

S 2 i = S? 2 = * r t* l5 (B.8c) 
where CJ has been absorbed into 3> r , <& r [7 — > <& r , and 



t = V^-P 2 - (B.9) 

Note that (|B.7[) implies that ||p|| < 1. 

Eq. (|B.8j) and (|B.9j) can be interpreted as follows: The component can be viewed 
as a discrete reflector sandwiched between two unitary transmission sections. The 
discrete reflector provides coupling between equal modes that propagate in opposite 
directions, and the unitary sections provide coupling between different modes in the 
same direction. For the discrete reflector, the reflection response from the left and 
right is p and — p, respectively, and the transmission is t. For the two unitary sections, 
there are no reflections, and the transmission responses from the left are 3>i and 
3? r , while the transmission responses from the right are and . Note that 
this interpretation is consistent with the reciprocity and lossless conditions (|B.4[) and 
(|B.5[) . for each of the three sections separately. By inspection, we find that (|B.8[) is 
invariant if PpP T — * p, PtP T — > t, P<&i — * 3>i, and <& r P T — > $ r where P is a 
real unitary matrix. Here P represents an arbitrary rotation of the eigenaxes of the 
reflector (p and t are now real and positive semidefinite). 

Using (|B.8p . the transfer matrix (|B.3[) can be written 



A* B* 
B A 

where the blocks A = *,*t~ 1 $;* and B = —&*t ; _1 p€>; satisfy 



(B.10) 



A* A — B T B* = I (B.lla) 
AB T — BA T = (B.llb) 
A T B* — B^A = 0. (B.llc) 
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